function mbs

deltheta = -45*pi/180;

% Try reducing the tolerance and observe the effect on the drift of total
% energy of the system
options = odeset('RelTol',1.0E-7','AbsTol',1.0E-7);

[t,x] = ode45(@dot,[0 10],[pi/2+deltheta 0]',options);

for i=1:length(t)
    [e(i),ke(i),pe(i)] = calcEnergy(x(i,:));    
end
hold off
plot(t,e);
hold on
plot(t,ke,'r');
plot(t,pe,'g');
legend('Total Energy','Kinetic Energy','Potential Energy');

figure(2)
plot(t,e-e(1));
legend('Drift in total energy is proportional to the integrator tolerance');

keyboard
function op = dot(~,x)

g = -9.81;
l = 1;
q = x(1,1);
qdot = x(2,1);

udot = -(g*cos(q)+cos(q)*qdot^2)/(2*l*(sin(q) + 1));

op = [qdot;udot];

function [e,ke,pe] = calcEnergy(x)

q = x(1,1);
qdot = x(1,2);

ke = (qdot^2*cos(q)^2)/2 + (qdot^2*(sin(q) + 1)^2)/2;
pe = -9.81*(sin(q)+1);
e = ke+pe;